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ABSTRACT 

We use the high resolution cosmological TV-body simulations from the Aquarius 
project to investigate in detail the mechanisms that determine the shape of Milky 
Way- type dark matter haloes. We find that, when measured at the instantaneous 
virial radius, the shape of individual haloes changes with time, evolving from a typ- 
ically prolate configuration at early stages to a more triaxial/oblate geometry at the 
present day. This evolution in halo shape correlates well with the distribution of the 
infalling material: prolate configurations arise when haloes are fed through narrow fila- 
ments, which characterizes the early epochs of halo assembly, whereas triaxial/oblate 
configurations result as the accretion turns more isotropic at later times. Interest- 
ingly, at redshift z = 0, clear imprints of the past history of each halo are recorded 
in their shapes at different radii, which also exhibit a variation from prolate in the 
inner regions to triaxial/oblate in the outskirts. Provided that the Aquarius haloes are 
fair representatives of Milky Way-like IO^^Mq objects, we conclude that the shape of 
such dark matter haloes is a complex, time-dependent property, with each radial shell 
retaining memory of the conditions at the time of collapse. 

Key words: galaxies: haloes, galaxies:formation, galaxies:evolution, cosmology:dark 
matter 



1 INTRODUCTION 

In our current understanding of the Universe, dark matter 
haloes constitute an integral part of galaxies. Their proper- 
ties, especially their density profile and shape, have received 
significant attention in recent years as they have been ar- 
gued to be sensitive to the fundamental properties of the 
dark matter particles. Numerical simulations have been ex- 
tensively used to study the characteristics of the dark matter 
haloes, exploring for example the effects of the environment, 
mass assembly history and the nature of dark matter itself 
(e.g. 



Bullock 2002 Bett et al. 2007 iMaccio et al. 2007 



|Hahn et al.|2007| [Spergel Steinhardti2000| [Yoshida et al!] 



* E-mail: cavera@astro.rug.nl 



2000) lAvila- Reese et al.l|200T] [Strigari et aL]|2007l [Wangfc] 



White||2009 ) 



The first fully analytical models of the formation of 
dark matter haloes such as the top-hat spherical collapse 



model (Gunn & Gott 



1972 



Fillmore & Goldreich 



1984), 



considered highly symmetric configurations. However, the 



pioneering work of Frenk et al. ( 1988 ); Dubinski & Carlberg 



(1991); Warren et al. (1992); Cole Lacey (1996); Thomas 



et al. ( 1998 ) demonstrated important deviations from spher- 



ical symmetry by measuring the shape of dark matter haloes 
in numerical A/'-body simulations evolved in a fully cosmo- 
logical context. These authors consistently found that after 
virialization, dark matter haloes are triaxial with more pro- 
late shapes towards the centre and more oblate shapes in 
the outskirts. Recent high resolution A^-body simulations 



2 Vera- Giro et al 



have yielded similar conclusions (|Jing & Suto||2002| Bailin 


& Steinmetz |2005| Kasun & Evrard 


|2005| IHopkins et al. 


2005 


Bett et aL"2007; Hayashi et a] 


. |2007( IKuhlen et al. 


2007 


Stadel et al. 2009; Diemand k Moore|2011|. 



Further studies based on numerical simulations have 
also revealed that the environment and mass assembly his- 
tory of a halo may play a crucial role in determining its 



shape. Pioneering work by |Tormen ( 1997) andjColberg et al. 
(1999 J suggested that the anisotropic infall of matter onto 
cluster-sized haloes was largely responsible for their shape, 
orientations and dynamics at different times. Because in- 
fall is governed by the surrounding large scale structure, 
we expect significant correlations between the halo shapes 
and their environment, although evidence both against and 
in support of such trends have been reported so far in the 
litera ture (|Lemson Ka uffmanrT 1999' "Avila- Reese et al.' 
[20051 IFaltenbacher et al.P2005 ; Altay et al. 2006 j Basilakos 
et al. 2006; Gottlober k Turchaninov 2006; Patiri et al.|2006| 
i Aragon-Calvo et al.,,2007, |Hahn et al.[(2007| [Maccio et al.| 
|2007| [Ragone-Figueroa et al.||2010| ). 

The observational determination of the shapes of dark 
matter haloes is challenging. Preferably dynamical tracers 
at large radii are to be used, but these tracers are by defini- 
tion rare. In external galaxies, constraints on these shapes 
have been put using the intrinsic shape of galactic discs 



(Fasano et al. 1993), the kinematics and morphology of the 



Hi layer ( 011ing^^l996j ^Becquaert Combes 1997^ ; ^Swaters^ 
et al. 1997), the morphology, temperature profile of X-ray 



isophotes fBuote & Canizares'1998( Buote et al.|2002 ), grav- 
itational lensing (Hoekstra et al. 2004) and the spatial distri- 



bution of galaxies within groups ( [Robotham et al . 2008 ) (for 
earlier reviews on the subject see | Rix||1996 Sackett^^l999J . 
The general consensus of all these studies is that haloes tend 
to be roughly oblate, with the smallest axis pointing perpen- 
dicular to the symmetry plane defined by the stellar compo- 
nent. Most of these constrains, however, pertain to the inner 
regions (a few optical radii) of galaxy-scale haloes. 

In the case of the Milky- Way, the shape constraints of- 
ten rely on the kinematics of individual stars, and include 
e.g. the use of the tilt of the velocity ellipsoid for nearby 
stars (|Siebert et al.||2008|, the proper motions of hyperve- 



locity stars ( Gnedin et al.||2005 ) or the dynamics of stellar 



streams ( Koposov et al.||2010 ). Interestingly, the use of the 
latter have provided contradictory results, notably in the 
case of the Sagittarius Stream. For example, the positional 
information was used to argue that the Milky Way halo is 
nearly spherical ( Ibata et al.||200l' Martmez-Delgado et al 



2004[ [Johnston et al.||2005| whereas the kinematics of stars 
in the leading stream could only be fit in a prolate halo 
elongated perpendicular to the disc ( Helmi|2004 Law et al 
2005|). More recently, [Law et al.] ( |2009| ); [Law &: Majewski 
(2010) have explored triaxial potentials with constant axis 



ratios, and found models that were able to fit simultaneously 
these constraints, albeit not completely satisfactorily. 

From the theoretical perspective, there is a gap in our 
knowledge about the way in which dark matter haloes ac- 
quire their shape, and in particular, on the impact of the 
dynamics of the surrounding large scale structure in the 
non- linear regime ( Lee et al.|2005 Betancort-Rijo & Trujillo 



2009| [Rossi et al.||2010| [Salvador-Sole et al.||2011| ). Indeed, 
most of the theoretical works cited above restrict their anal- 
ysis to the present day correlations with the environment. 



and do not consider when the shapes have been established 
and how they relate to the past history of an object. In this 
paper we use state-of-the-art high resolution A^-body sim- 
ulations that track the formation of five ^ lO^^M© Milky 
Way-like haloes in a fully cosmological context, in order to 
gain further understanding of this problem. 

This paper is organized as follows. In Section [2] we de- 
scribe the simulations used in this work which are part of the 
Virgo Consortium's Aquarius project, and test in Section |3] 
the convergence of halo shapes for different resolutions. The 
shapes at the present day, evolution and their relation with 
the formation history of the dark matter haloes are analyzed 
in Section |4] and [5] Finally, we discuss and summarize our 
main conclusions in Section[6] For completeness, we compare 
in Appendix [a] the results obtained by using several differ- 
ent schemes to determine halo shapes. Also a novel method 
for measuring the size of the filaments based on dynamical 
arguments is described in Appendix [B] 



2 NUMERICAL SIMULATIONS 

We use the Aquarius Simulations, a suite of high resolution 
A/'-body cosmological simulations of six Milky Way sized 
dark matter haloes (Springel et al.|2008 ). These haloes were 
selected from a larger ACDM cosmological box of 100h~^ 
Mpc side with parameters Qrn = 0.25, Qa = 0.75, ag = 0.9, 
= 1 and Ho = lOO/i km s"^Mpc"^ = 73 km s"^Mpc"\ 
The Aquarius haloes, labeled Aq-A to Aq-F, have a final 
mass ^ lO^^Mo and were chosen to be relatively isolated 
at redshift z = 0. The selection procedure was otherwise 
random, which allows us to study the impact of different 
assembly histories on the shape of Milky Way sized dark 
matter haloes. 



Each halo has been re-simulated at various resolution 
levels that accurately replicate the power-spectrum and 
phases for the resolved structures in all runs. Following the 
notation introduced in previous papers, we refer to each level 
of resolution as -5 to -1, for the lowest to highest resolution. 
The mass per particle varies from rup = 2.94 x lO^/i"^M0 
in the level 5 to rrip = 1.25 x lO^/i^^M© for the highest res- 
olution run, which resolves a given halo with approximately 
half-million up to 1.5-billion particles within the virial ra- 
dius for level 5 and 1, respectivel}Q The results discussed 
in this paper pertain mostly to the 4th level of resolution 
(nip ^ 2 X lO^/i~^M0, rivir 10^ particles within rvir, grav- 
itational softening e = 250h~^ pc) of haloes Aq-A to Aq-E 
and aim to target the possible structure of the Milky Way 
dark matter halo. The Aquarius halo Aq-F is not considered 
in this work because it experiences a major merger less than 
~ 5 Gyr ago; this is unlikely to be consistent with our cur- 
rent view of the assembly of the Milky Way halo, expected 
to have had no major mergers since z ^ 1 (Toth & Ostriker 



1992) 



In this context, this paper focuses on dark matter haloes 



^ Throughout this paper we refer to the virial radius, rvir, as the 
spherical radius that contains a mean density equal to 200 times 
the critical density of the Universe at a given time. Other virial 
quantities, such as mass and velocities, mvir and i;vir respectively, 
correspond to those measured within rvir- 
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Halo 


nip 


Tvir 


TTT-vir 


nvir/lO^ 


''"conv 


Aq-A-4 


2.87 


179.36 


1.34 


4.68 


2.27 


Aq-B-4 


1.64 


137.86 


0.61 


3.72 


2.14 


Aq-C-4 


2.35 


177.89 


1.31 


5.58 


2.07 


Aq-D-4 


1.95 


177.83 


1.31 


6.69 


2.14 


Aq-E-4 


1.90 


155.95 


0.88 


4.64 


2.05 



Table 1. Numerical details of the Aquarius haloes for resolution 
level-4. Each column lists: (1) the name of the halo, (2) mass per 
particle rup in 10^ , (3) the virial radius rvir in kpc 
and (4) virial mass mvir in IO^^Mq at 2; = 0, (5) rivir the 

number of dark matter particles within rvir, (6) rconv the con- 
vergence radius in kpc h~^. The gravitational softening length is 
e = 2b0h~^ pc for all haloes at this resolution level. 




10.0 
R (kpc h"^) 



100.0 



Figure 1. Convergence in the axis ratios for the Aq-A halo at 
four different numerical resolutions. The agreement is remarkable 
showing that the method used provides reliable axis ratios down 
to the "convergence radius" rconv (denoted by the vertical lines). 



with a relatively quiescent merger history (mass ratios lower 
than 1:5) after z ^ 2 (for a detailed comparison of the 
Aquarius haloes' mass accretion history with respect to the 
general expectations for haloes of similar virial mass see 
Boylan-Kolchin et al. 2010). We summarize the most rel- 
evant numerical details of our haloes at level 4 in Table 111 



and we refer the interested reader to Springel et aL| (|2008| 
for further information. 



3 HALO SHAPE DETERMINATION AND 
CONVERGENCE 

Various methods have been introduced in the literature to 
measure the shapes of dark matter haloes. In Appendix [A] 
we describe a wide range of these methods in detail, and 
show that the measured halo shapes agree reasonably well 
when applied to the same object. Throughout this article 
we use the "reduced" inertia tensor method as implemented 
recently by [Allgood et al.| ([2006| . This tensor is defined as 




Aq-A-4 



500h-' kpc 




Figure 2. Present day dark matter distribution of the Aquarius 
haloes and their surroundings. Colors are proportional to the log- 
arithm of the squared dark matter density integrated along the 
line of sight. In each panel the depth of the projected image is 
1.5h~^ Mpc. The virial ellipsoid is indicated in each panel with 
a solid white line and the green arrow shows the direction of the 
minor axis of the halo at the virial contour. The system is orien- 
tated such that the minor axis points horizontally and the major 
axis points vertically. 



E 



X 1 X 



dl 



(1) 



where dk is a distance measure to the k-th particle and V 
is the set of particles of interest. Assuming that dark mat- 
ter haloes can be represented by ellipsoids of axis lengths 
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a ^ b ^ c, the axis ratios q — h/a and s — c/a are the 
ratios of the square-roots of the eigenvalues of /, and the 
directions of the principal axes are given by the correspond- 
ing eigenvectors. Initially the set V is given by all parti- 
cles located inside a sphere which is re-shaped iteratively 
using the eigenvalues of /. The distance measure used is 
d\ — x\-\-y^/q^ -\- z^/ ^ where q and s are updated in each 
iteration. Furthermore, in the figures below we have removed 
all bound substructures contained in a halo using the Sub- 



find algorithm ( [Springel et al. 2001). This alleviates the 
noise and artificial tilting of the ellipsoids that is introduced 
by such substructures. More details on our implementation 
can be found in the Appendix [A] 

We test the convergence of the shape measurements us- 
ing the different resolutions of the halo Aq-A, from level 5 to 
level 2. There is an increase of a factor ^ 230 in the number 
of particles within the virial radius and a factor ~ 10 reduc- 
tion in gravitational softening from Aq-A-5 to Aq-A-2. Fig.^ 
shows the axis ratios as a function of R for halo Aq-A, where 
R is defined as the geometrical mean of the axis lengths 
R — {ahcY^^ . This quantity then provides a notion of dis- 
tance to the centre of the ellipsoid with the advantage of 
being volume invariant. Hereafter, we use R when referring 
to quantities measured using ellipsoidal contours, and r for 
spherical contours. With this convention, rvir is the spheri- 
cal radius enclosing an average overdensity of 200 times the 
critical value, whilst i?vir refers to the ellipsoid enclosing the 
same overdensity. 

The agreement between the different resolutions in 
Fig. [l] is remarkable at all radii, although small deviations 
are present in the inner regions where resolution effects are 
expected to be important. The vertical lines in the lower 
panel indicate the "convergence radius" rconv (Power et al. 



2003 Navarro et al. 2010| ). This radius is defined by set- 
ting /^(rconv) = trelax(r'conv)/tcirc(rvir) = 7, whcrC trelax is 

the local relaxation time and tcirc(^vir) corresponds to the 
circular orbit timescale at rvir- This choice ensures that cir- 
cular velocity profiles for r > rconv deviate less than 2.5% 
from the highest resolution value. Fig. [l] shows that rconv 
provides also a very good estimate for the smallest radius 
at which the halo shape has converged. Although not ex- 
plicitly shown, we have tested that the orientation of the 
ellipsoids is a well defined attribute independent of resolu- 
tion for rconv < R < Rvir- 

Our analysis suggests that convergence in halo shapes 
is achieved in the level 4 from approximately 2h~^ kpc on- 
wards (Table [T]). We will therefore focus in what follows on 
the study of the different Aquarius haloes Aq-A to Aq-E at 
the level of resolution 4, since this yields a robust charac- 
terization of the dark matter halo shapes and orientation, 
while keeping the numerical cost constrained. 



4 THE SHAPE OF DARK MATTER HALOES 
4.1 Present day 

Fig. [2] shows a snapshot view of the Aquarius haloes and 
their surrounding environment. The color scale used is pro- 
portional to the logarithm of the squared dark matter den- 
sity integrated along the line of sight, with a projection 
depth per panel of 1.5/i~^Mpc. All images have been rotated 



according to the orientation of the virial contours (indicated 
by the solid white curve) in such a way that the major and 
minor axis define the vertical and horizontal direction, re- 
spectively. 

A careful look at the environment of each object re- 
veals that the Aquarius haloes are all embedded within a 
filamentary-like structure, which is better defined in some 
cases than in others (e.g. Fig. [2] Aq-A vs Aq-E), a result that 
holds independently of the projection used. These filaments 
also seem to contain the most massive substructures in the 
box and interestingly, the minor axes of the virial contours 
tend to be perpendicular to the direction defined by these. 
This is in very good agreement with the statistical find- 



ings of Bailin & Steinmetz (2005), who analyzed a sample 
of ^ 4000 haloes in a wide range of masses. One relevant 
consequence of this type of configuration is that most of the 
substructure will preferentially be accreted along the major 
axis of the halo, a feature that has been used to explain 
the preferential alignment of satellite galaxies with respect 
to their hosts in observational studies of external galaxies 
(e.g. Brainerd""2005', and references therein), in groups (e.g. 
Binggeh 1982; Kitzbichler Saurer 2003 i Yang et al. 2006, 



2011,) and also numerical simulations (e.g.|Knebe et al.|2004 


Libeskind et al.| 


2005, 2007; Kang et al.||2007| |Sales et al. 


20071 IBailin et a 


1I2OO8I IFaltenbacher et al.|2008| |Libeskind 


et al.||2011|). 



The halo shapes and orientations as a function of radius 
are shown in Fig. [3] On the right panel, besides the axis 
ratios, we also plot the triaxiality parameter T defined as 



(2) 



which is unity for a perfect prolate distribution and zero in 
the oblate case ( |Franx et al. 1 1991 1 [Warren et al.|1992t . This 
figure shows that in the inner regions all dark matter haloes 
are more "prolate" with b/a ^ c/a ^ 0.4-0.6 ( [Hayashi et al.| 
|2007| ). On the other hand, at large radii b/a increases to 0.8- 
0.9, and the mass distribution becomes more oblate/triaxial. 
At intermediate radii haloes are typically triaxial, but no- 
tice that the radius at which the prolate-to-oblate transition 
occurs (defined as the radius where T = 0.5) is different for 
each object, ranging from ^ 20h~^ kpc for Aq-E to ^ 100/i~^ 
kpc in the case of Aq-A. 

Regardless of the change in shape of the concentric ellip- 
soids, these tend to remain well aligned throughout the halo. 
The right panel of Fig.|3]shows the cosines of the relative an- 
gle between the major, intermediate and minor axis of the 
ellipsoids at different radii and at the virial contour. The 
alignment is remarkable, with the exception of halo Aq-D, 
where the intermediate and minor axis in the inner regions 
are rotated ^ 60° with respect to their orientation at the 
virial radius, an issue that we explore in the next section. 



4.2 Evolution 

Fig. ^ shows a series of snapshots of halo Aq-D at different 
stages of its evolution. The orientation of the coordinate sys- 
tem is now kept fixed for all the snapshots. The correspond- 
ing time is quoted in the upper left corner of each panel, the 
virial ellipsoids are depicted by white ellipses while the pro- 
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Figure 3. Axis ratios (left) and directions (right) as a function of K for each of the Aq-haloes. In the left panels the thick solid lines 
represent the axis ratios hja (red) and c/a (blue), while the dashed curves are the triaxiality parameter T. In general haloes are more 
prolate in the inner parts and more oblate (and triaxial) in the outskirts. All concentric shells are also strongly aligned with the exception 
of halo Aq-D which exhibits a clear twisting close to the virial radius. 



jection of the minor axes onto the plane are indicated by the 
green arrows. There are several points worth highlighting: 

• The shape and orientation of the halo seem to change 
throughout time (with the caveat that we are just seeing its 
evolution in projection). 

• Also evident from this figure is the filamentary struc- 
ture that characterizes the surroundings of Aq-D. Notice that 
the filament where this halo is located at the present day 
was already in place at t < 1 Gyr, and fully dominates the 
environment from t ^ 3 Gyr onwards, with remarkable co- 
herence in time and direction. 

• The relative orientation of halo Aq-D with respect to 
its environment is interesting since its minor axis is parallel 
to the direction of the filament at late times. This seems 
to be at odds with expectations from statistical studies of 



dark matter halo alignments (Bailin & Steinmetz 20051 



Fal- 



tenbacher et al.|2005; Aragon-Calvo et al. 2007; Zhang et al. 



2009| ) as well as in contradiction with the other haloes illus- 
trated in Fig. [2] and the analysis presented in Section [XT] A 
closer inspection of the history of this halo shows that the 
infall of material at t > 7.5 Gyr occurs mostly along a sec- 



ondary filament (only barely visible in this projection) whose 
direction is almost perpendicular to that of the most promi- 
nent and massive filament that dominates the surrounding 
large scale structure and is clearly seen in Fig. [4] This change 
of infall direction explains the apparent change in orienta- 
tion of halo Aq-D at late times. 

• The relative size of a filament with respect to that of 
the halo increases with time: the filament's cross section 
is comparable to the virial radius at t < 5 Gyr, whereas 
the dark matter halo becomes smaller than and is fully em- 
bedded in the filament at later times. This has interesting 
consequences on the halo shape as we discuss further in Sec- 
tion [5] We have explicitly checked that this conclusion is not 
dominated by projection effects nor rendering of our images. 
We refer the reader to Appendix ^ where we introduce a 
physically motivated definition of a filament (based on the 
number of caustic-crossings a particle has experienced) and 
compare its size with that of the halo at each timestep. 

In Fig. [5] we quantify the time evolution of the axis ra- 
tios measured at the instantaneous virial ellipsoids for haloes 
Aq-A to Aq-E. For reference, we have indicated with grey la- 
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Figure 4. Time sequence of the formation of halo Aq-D from t ^ \ Gyr to the present day. Each box shows the dark matter density 
distribution around this halo together with the shape and orientation of the virial contour at the given time. As before, the green arrows 
indicate the projection of the minor axis onto the plane, with the length of the arrow being proportional to the minor axis length. The 
shape of this halo evolves continuously with time, showing stretching and tilting in response to the infall of matter. The halo is oriented 
along the filament up to t ~ 7 Gyr and then rotates and becomes perpendicular to it at the present day. This figure clearly illustrates 
the complexity involved in the build up of the present day dark matter halo's structure. 



bels the corresponding physical sizes at a given time. We will 
only follow the evolution from i — 2 Gyr onwards (which 
corresponds roughly to the time at which the haloes have 
~ 20% of their final mass), which guarantees that the cen- 
tre of mass as well as the shapes at the virial radius are well 
defined at all times for the full sample of haloes. 

Fig.|5]shows that the shape of the virial ellipsoids is not 
constant with time. In general, haloes seem to evolve from a 
quite prolate configuration at early times towards more tri- 
axial shapes at the present day. In most of our haloes (with 
the exception of Aq-D), this evolution in shape is driven by 
a larger increase of the intermediate-to-major axis ratio hja 
than in c/a. Except for this weak general trend, the evolu- 
tion of the axis ratios is quite disparate from halo to halo. 

A qualitative comparison between the left panel of 
Fig. [3] and [5] suggests a certain degree of resemblance be- 
tween the overall evolution of the virial ellipsoid's shape with 
time (as measured, for example, by the triaxiality parame- 
ter) and the shapes measured as a function of radius (Fig.[3| 
for each of our haloes, despite their intrinsic differences. This 
suggests that the dark matter haloes retain certain memory 
of their configuration in the past, and that this is imprinted 



on their present day structure. We note however that this 
refers to the overall trend with radius/time and does not im- 
ply that it is possible to recover the exact numerical value of 
the axis ratios at a given time from their present day value 
at a given radius. 

This analogy between radius/time is better seen in 
Fig.[6] where we plot the present day axis ratios hja vs c/a at 
different radii (grey curve) as well as those measured at the 
virial ellipsoids for different times (color points) . Because we 
sample from t — 2 Gyr onwards, a dotted or a solid line is 
used to indicate respectively, the shapes at radii smaller or 
larger than the virial ellipsoid i?vir at this time. This figure 
shows that haloes evolve away from the 1:1 line and there- 
fore tend towards less prolate shapes with time (and also 
radius), with the clear exception of halo Aq-D. Moreover, 
the reasonable agreement between the solid curve and the 
diamonds in this figure suggests that the present day shape 
of the dark matter haloes at different radii provides infor- 
mation about the evolution of the virial ellipticities during 
their assembly history. 

The analysis presented in this section, if generalized to 
Milky Way- type haloes, is directly relevant to the modeling 
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Figure 5. Temporal evolution of the axis ratios and triaxiality 
parameter for the virial contour. The color-coding is the same as 
in Fig. [3] A comparison to that figure shows that a correlation 
exists between the shape as a function of time, and the present- 
day shape as a function of radius. To aid this comparison we have 
added a second axis (gray), indicating the size of the virial contour 
i^vir- This is derived for each halo by fitting: -Rvir(^) = A{1 — 
e(*-'^)/(2*h)) to the measured evolution of the virial (ellipsoidal) 
radius. Here is the time at which the virial mass reaches half 
its present value, and A and r are free parameters. 



and interpretation of observational constraints on the shape 
of the Galactic potential. As discussed in the Introduction, 
the Sagittarius stream has led to contradictory results when 
modeled in axisymmetric dark matter haloes with constant 
axis ratios (in t he potential) (|Ibata et al.|2001| |H elmi"20Q4"; 
Law et al.|2QQ5| [Johnston et al.|2005||Fellhaue"ret al . 2006). 
We have shown here that significant variations in the axis 
ratios as a function of radius exist for all our five Aquar- 
ius haloes and those variations are linked to the evolution- 
ary history of each object. Although this might add an ex- 
tra degree of freedom to models that attempt to constrain 
the Galactic potential, our results from Fig. [6] indicate that 
by doing so it may be possible to retrieve the local condi- 
tions around the Milky Way's halo throughout its assembly 
( Banerjee k Jog||2011 ). 



5 ENVIRONMENT, MASS ACCRETION AND 
HALO SHAPES 

Dark matter haloes continue to grow with time, through a 
regular, although not always steady, injection of matter. Ha- 
los will respond to this new material and re-structure them- 
selves onto a new equilibrium configuration as part of the 
virialization process. The environment of an object deter- 
mines to some extent the way in which the material is incor- 
porated into the halo. For instance, thin filaments imply ac- 
cretion through well-defined preferential directions whereas 
a more isotropic mode is expected when the halo is embed- 
ded in a larger structure. Several correlations of halo shapes 
with environment have been found so far in simulations (e.g. 



Avila- Reese et al.'2005';'Faltenbacher et al.'2005';'Patiri e tliL 
2006: Basilakos et al. 2006; Maccio et al. 2007; Hahn et~ 



2007 



Flin 



and observations (e.g. [Pimbblet |2005| [Godlowski k, 
2010| [Niederste-Ostholt et al.||2010| |Paz et al.||20lT^ 



albeit with significant scatter. In this section we explore in 
detail the role played by the environment and mass assembly 
histories in shaping the Aquarius haloes. 

The angular distribution on the sky of the infalling mat- 
ter can provide useful information about the preferred di- 
rections and modes of the ongoing accret ion (Tormen|1997| 
Colberg et al.||1999| |Aubert et al.||2004| |Aubert Pichoii| 
Libeskind et al.|20TT ). For instance, whereas isotropic 



2007 



accretion would lead to a uniform signal on the sky, the 
presence of a thin filament will give rise to a bi-modal dis- 
tribution of points in two opposite (180° apart) directions. 
This can be quantified further by means of a multipole ex- 
pansion of the infalling particles on the sky at a given time. 
We measure the power spectrum for the mode / as. 



Ci 



1 1 \^rn 



(3) 



m— — l 

where the expansion coefficients are 



(4) 



with Qfc being the angular position of the k-th particle cross- 
ing the virial radius rvir with negative radial velocity at any 
given time. The asterisk on the right hand side of Eq. Q 
indicates the complex conjugate of the term within paren- 
thesis. 

In the scheme introduced above, the / = term 
(monopole) is a constant equal to unity and used for the 
overall normalization of the expansion. Notice that al- 
though this choice is arbitrary, the relative relevance of the 
monopole with respect to all the other modes, Co/T^iCi, is 
an indication of how isotropic the distribution is; i.e. a per- 
fectly isotropic infall would correspond to all the power in 
the / = term. On the other hand, a significant contri- 
bution of the / = 2 (or quadrupolar moment) term arises 
when the accretion occurs through a well-defined direction 
in space, i.e. a filament. Similarly, accretion corresponding 
to more than one preferential direction will shift the power 
away from / = 2 and towards higher moments. Notice that a 
point mass on the sky will excite, by definition, a wide range 
of modes with similar power; just like the Fourier transform 



8 Vera- Giro et al 




Figure 6. Axis ratios as a function of position (gray line) and time (diamonds). The colors indicate the time when the axis ratios at 
the virial contour have been measured. Notice that since we consider temporal evolution only after t = 2 Gyr, we have differentiated the 
shape as a function of position (gray curve) using solid (dotted) lines for distances larger (smaller) than the virial radius of each halo at 
t = 2 Gyr (see text for details). 



of a Dirac-delta function has constant power for all modes in 
the frequency space. We therefore expect single satellite in- 
fall events to excite higher moments compared to the smooth 
accretion. In the limit of a satellite that occupies a large area 
of the sky, the configuration will then resemble a dipole and 
the power spectrum will exhibit high power in the / = 1 
mode. 

In practice, most of the information is encoded in the 
terms / ^ 2 (e.g. |Quinn Binney|1992 [Eisenstein Loeb] 



1995 Aubert et al. 2004"). We have checked in our exper- 



iments that higher moments than quadrupolar contribute 
always less than ^ 15% of all power at any given time once 
substructures have been removed. We therefore focus our 
analysis mostly on the / = and / = 2 moments since they 
provide most of the information that allow to characterize 
the mass accretion onto dark matter haloes. 

Figj7| shows this multipole decomposition introduced 
in Eq. ([3| of the mass infalling onto Aq-A-4 as a function of 
time. At each output time, we select particles with negative 
radial velocity (infalling), Vr < 0, that are in the spherical 
shell 1.0 ^ r/rvir ^ 1-2 and compute the corresponding Ci 
of the distributiorj^ The upper and the middle panels of 
this figure show the Ci power spectra obtained respectively, 



by including and by removing particles associated with sub- 
structures as identified by Subfind. Both distributions are 
quite similar, although as expected from the discussion in 
the previous paragraph, satellite accretion excites in general 
higher modes that last a very short timescale, and are visible 
as clear "spikes" in the upper panel of the figure. 

The virialized regions of galaxy-sized objects can ex- 
tend well beyond their formal rvir (e.g. Cuesta et al.|2008 ), 
introducing a signal in the power spectrum that is driven 
by the halo's intrinsic shape rather than by the surrounding 
infall pattern. In order to avoid confusion with the mate- 
rial already in place and in equilibrium in the outskirts of 
the halo, we selected also the subsample of particles within 
1.0 ^ r/rvir ^1.2 that are infalling for the first time onto 
the halo (in practice we do this by requiring that a particle 
in a given output has never been within the virial radius of 
the main object at any previous time step). The correspond- 
ing power spectrum (after removing the subhaloes' contri- 
bution) is shown in the bottom panel of Fig. [t] This distri- 
bution agrees well with those shown in the other panels of 
this figure, although the features appear noisier due to the 
smaller number of particles being considered. 

The left panel of Fig. [S] shows the power spectrum of 



^ We have tested that the qualitative behavior shown in this fig- 
ure and the relative relevance of each mode with respect to the 



whole spectrum does not depend on the particular shell that is 
analyzed for 1 < r/rvir < 2 for all of our haloes. 
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all infalling material for all five Aquarius haloes after the 
contribution from subhaloes has been removed. The resid- 
ual "spikes" visible in this figure correspond to the matter 
that is associated to infalling substructures but that is rather 
loose and, consequently, has not been assigned to a partic- 
ular subhalo by the halo finder. 

The right panel of Fig.[8]shows the relative contribution 
of the / = 2 (solid black curve) and / = (blue dot-dashed) 
moments to the total power spectrum as a function of timaj 
These have been computed using only the subset of parti- 
cles on their first infall (although these curves do not change 
significantly when all infalling particles are considered in- 
stead). Large C2 values are associated with the presence of 
net filamentary accretion. Fig.[8]shows that this condition is 
typically found at early times in our haloes, with the excep- 
tion of Aq-B, which shows no clear sign of smooth accretion 
through a filament at any time. The halo Aq-D shows also 
a peculiarly high power in the / = 2 mode at late times, 
whereas for most of the objects the relevance of C2 remains 
approximately flat (and negligible) in the last few Gyr (a 
feature which is also evident in the left panel of Fig. [8|. 

The monopole term largely dominates the accretion at 
later times (blue dot-dashed curve in Fig. [sj in all haloes. 
Exceptions are haloes Aq-A and Aq-D, which show a decline 
in Cq/T^Ci beyond t 12 Gyr and t ^ 1^ Gyr respec- 
tively. In the case of Aq-A this is related to an increase in 
power of / ^ 2 modes, while for Aq-D this is driven only 
by the / = 2 moment, as mentioned in the previous para- 
graph. As we show in Appendix [B] this more isotropic infall 
is a direct consequence of the increase in the relative size of 
the filaments feeding the dark matter haloes: infalling par- 
ticles cover wider angles on the sky which leads to larger Co 
contributions with negligible / = 2 component. Notice that 
halo Aq-E is, in some sense, the most extreme object, with a 
monopolar term that dominates the power spectrum of infall 
material during almost all its entire history {Cq/T^iCi > 0.9 
for t ~ 3 Gyr onwards). 

A careful comparison of Fig.[8]and Fig. [5] reveals a good 
correlation between the infall of material through filaments 
(large C2) and the shape of the haloes at a given time: high 
/ = 2 moments are associated with virial contours that turn 
prolate (e.g. Aq-A for t < 6 Gyr, Aq-C for t < 4 Gyr, Aq-D 
for t > 8 Gyr). In particular, this multipole decomposition 
shows a filamentary accretion mode for halo Aq-D that is 
present at t ^ 10 Gyr, helping to explain why its axis ratio 
b/c remains close to unity at quite late times. When the 
accretion is more isotropic (i.e. Co/^iCi :^ 0.9 ), the haloes 
become most nearly oblate. This explains why halo Aq-E is 
the more oblate in our sample, with b/a ^ 0.9 from t ^ 3.5 
Gyr onwards, and also why it has this shape at smaller radii, 
compared for example to halo Aq-C. 

We found that the injection of material that occurs 
along filaments leads to a more prolate halo shape, an effect 
that is naturally enhanced at early times due to the relatively 
smaller size of the filaments with respect to the dark mat- 



^ Note that we are not explicitly showing the I = term in Fig. [s] 
because it is, by definition, set to unity in our formalism. However, 
the "relative" importance of the monopole with respect to the 
contribution of all other moments is a well defined quantity, that 
we analyze in more detail on the right panel of the same figure. 



u" -4 3 -3.3 -2.3 ."^-^ 




Figure 7. Multipole expansion of the infalling material (vr < 
0) in the region 1.0 ^ r/rvir ^ 1-2 as a function of time for 
halo Aq-A-4. The upper and middle panels show the results with 
and without the contribution from subhaloes, respectively. In the 
bottom panel only those particles that are on their first infall 
have been considered. 



ter halo ( Avila- Reese et aL]|2005 Gottlober Turchan inovj 
2006 ). On the other hand, the more isotropic mass accretion 
that characterizes later phases of halo assembly for lO^^M© 
Milky Way-like objects, yields more oblate/triaxial geome- 
tries. When combined with the memory effect alluded to 
in Section [4. 2 1 we find that prolate shapes are naturally ex- 
pected to be set in the earliest collapsed regions (small radii) 
of ~ lO^^M© haloes, in good agreement with previous work 
dBailin Steinmetz|2005l|Hayashi et al.|2007t . The filamen- 
tary structure typically found at early times thus seems to 
be responsible for the general trend of our Milky Way-like 
haloes to be more prolate in their inner regions. 



6 CONCLUSIONS 

In this paper we analyzed the shape of five Milky Way-like 
dark matter haloes selected from the Aquarius A/'-body sim- 
ulations. We compared the performance of several methods 
proposed in the literature to measure the shapes of haloes 
and found good agreement between all techniques, especially 
in the inner regions where substructures play only a minor 
role. Using an implementation of the normalized inertia ten- 
sor algorithm described in Allgood et al. (2006), we have 
found excellent convergence between the several resolution 
levels of Aquarius, where the shapes can be robustly mea- 
sured down to the convergence radius. 

We find that mass assembly and environment are both 
responsible for setting the shapes of dark matter haloes. The 
early evolutionary phases of lO^^M© Milky Way-like haloes 
are characterized by the accretion of matter through narrow 
filaments. In these circumstances the haloes -as measured 
by their virial contours- are prolate and their minor axes 
tend to point perpendicular to the infall (filament) direction. 
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Figure 8. Left: Multipole expansion of the infalling material in the region 1.0 ^ r/rvir ^ 1-2 as a function of time once the contribution 
from subhaloes has been removed for the five Aquarius haloes. Right: Relative contribution of the / = (blue dot-dashed) and I = 2 
(black solid) modes to the total power spectrum as a function of time for particles on their first infall onto each halo. While 
provides information about the material infalling along a filament, large contributions from the monopolar term, Cq/TiCi, imply that the 
accretion is isotropic. Notice the good correlation between the time intervals with a clear signature of mass accretion through filaments 
(high power in ^ = 2 mode) and the more prolate shape in Fig. [s] 



Nonetheless, temporary tilting of the virial ellipsoids may 
occur when mass is accreted from a different direction. The 
latter is the case for just one of our haloes located in a well 
defined filament at redshift z — 0. 

On the other hand, at later times the cross-section of the 
filaments becomes larger than the typical size of Milky Way- 
like haloes and as a result, accretion turns more isotropic 
and the objects evolve into a more oblate/triaxial configu- 
ration. This transition does not occur at the same time for 
all haloes in the explored mass range but is strongly deter- 
mined by their individual history of mass assembly and their 
surrounding environment. 

The geometrical properties of haloes at different epochs 
are not lost: haloes retain memory of their structure at ear- 
lier times. This memory is imprinted in their present-day 
shape trends with radius, which change from typically pro- 
late in the inner (earlier collapsed) regions to a triaxial in 
the outskirts (corresponding to the shells that have collapsed 
last and are now at the virial radius). These results are in 



excellent agreement with previous findings (Bailin & Stein- 
metz||2005| [Hayashi et al.|2007t . 



A corollary of our results is that the strong link between 
halo properties and assembly history, which can show large 
variations from halo to halo, will make any instantaneous 
correlation between haloes shape or orientation and mass or 
environment rather weak, explaining in part the relatively 
large scatter in such trends found in earlier studies (e.g. 
Bailin Steinmetz||20Q5| |Bett et al.||2Q07 ). 



It is important to stress that we have neglected the 
effects induced by the presence of baryons on the dark mat- 



ter halo shapes (|Kazantzidis et al.||2004 Bailin et al.||2005 


Gustafsson et al. '2006 i 


Debattista et al.|2008||Pedrosa et al.| 


2009: |Lau et al. 2011 


Valluri et al.||2010|. Therefore our 



results may only be directly applicable to dark matter dom- 
inated objects such as low surface brightness galaxies. How- 

( [20T0t ;|Abadi et al.|([20T0) 



ever, the work of Tissera et al 



jBett et al. ( 2010 ) shows that even though the halo axis ratios 
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increase when a disc is formed (i.e. they become rounder), 
the trends with radius appear to be preserved. 

With the caveat of the neglected baryonic effects and 
the relatively low number of objects studied, our findings 
may be directly relevant to the modeling of stellar streams 
used to determine the gravitational potential of the Milky 
Way. For example, the inconsistencies found between the 
constraints imposed by the positions and by the kinematics 
of stars in the Sagittarius stream could be indicative of a 
change in the shape of the Milky Way halo with radius (al- 
though seelLaw et al.||2009| iLaw & Majewski||2010|. Such 



a change could in principle be measured by using stellar 
streams which are on different orbits in the Galactic halo. 
Furthermore, it may even be possible to employ these to 
determine the growth history and early environment of the 
Milky Way. A validation of these ideas, however, is bound 
to a proper evaluation of the effect of the baryonic matter 
on our findings, as well as to the study of larger statistical 
samples, both issues that we plan to address in the near 
future. 
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APPENDIX A: THE SHAPE OF DARK 
MATTER HALOES: METHODS 

In this section we provide a review of the various methods 
that have been proposed in the hterature to measure dark 
matter halo shapes and compare the results when applied to 
the same object. The methods most commonly used are the 
diagonalization of the inertia tensor and the characterization 
with ellipsoids of either the interpolated density field or the 
underlying gravitational potential (Warnick et al. 2008). We 
introduce as well an additional scheme, which incorporates 
the advantages of different pre-existing methods, and where 
halo shapes are determined by fitting ellipsoids to the 3D 
iso-density surfaces. 

Al Inertia tensor 

One of the drawbacks of methods based on the determina- 
tion of the inertia tensor is its quadratic relation with dis- 
tance, which assigns the largest weight to particles residing 
far away from the centre, an issue complicated further by 
the presence of substructures in the outskirts of the dark 
matter haloes. 

The bias introduced by this distance weighting can be 
alleviated by normalizing each coordinate with some mea- 
sure of distance ([Gerhard 1983( |Dubinski Carlberg|1991 ). 
In this case the "reduced" inertia tensor 



(Al) 



where dk is a distance measure to the k-ih particle and V is a, 
set of particles' positions. Assuming that dark matter haloes 
can be represented by ellipsoids, the axis ratios are the ratios 
of the square-roots of the eigenvalues of /, and the directions 
of the principal axes are given by the corresponding eigen- 
vectors. To determine the axis lengths (e.g. b and c) however 
requires knowledge of the third axis (a). Therefore there are 
different choices to be made in this method, namely, which 
is the initial set of points V, the distance measure d and the 
way in which a is defined. Different approaches have been 
followed to set these quantities. 



Warren et ah] ( |1992D use an iterative scheme keeping 



the ellipsoid volume constant. The set V, initially chosen to 
be a spherical shell, is iteratively deformed and reoriented 
using the eigenvalues and eigenvectors of the reduced inertia 
tensor. These authors take d to be the Euclidean distance to 



given particle. Bailin & Steinmetz (2005), however, argue 



that this systematically overpredicts the roundness of the 
haloes. In the figures below we have modified Warren et al.[ s 
method by using d^ = x\ + y1 / + zl / instead, and where 
q — hja and s — cja are updated in each iteration. We define 
also the shell's radii using this definition of distance, instead 
of the Euclidean one. 
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Figure A2. Axis ratios (a ^ 6 ^ c) for Aq-A-4 as a function 
of the ellipsoidal distance R = (abcy^^. The solid red line cor- 
responds to our new algorithm, while dot-dashed black, dashed 
cyan, dot-dashed blue and dotted green correspond respectively, 
to: iso-potential conto urs ([Hayashi et al. 2007| >, particle-based iso- 
density inertia tensor ( |Jing Suto|2002 |, normalized inertia ten- 
sor diagonalized upon hollow ( [Warren et al.| |l992) and on solid 
( [Allgood et~ar|[2006[ > ellipsoids. Vertical arrows indicate the size 
of the virial ellipsoid, R^ir, for each method. All methods agree 
well, especially in the inner regions (R < 100h~^ kpc). 



Allgood et al. (2006) also employ an iterative scheme 



but now keeping the largest axis length constant. Initially 
the set V is selected to be given by all particles located inside 
a sphere (as opposed to a spherical shell of a given radius) 
which is reshaped iteratively using the eigenvalues. As be- 
fore, the orientation is determined from the eigenvectors of 
/, and the distance measure used is dl = xl -\- yl / q'^ -\- zl / . 
In the figures below we have removed all bound substruc- 
tures contained in a halo. This alleviates the noise and arti- 
ficial tilting of the ellipsoids that is introduced by such sub- 
structures. We use SUBFIND j Springel et al.|200T ) to identify 
and remove particles associated to subhaloes within the re- 
gion of interest. 

In our implementations the iterations are stopped when 
convergence in the axis ratios is reached, which we take to be 
when the relative change is smaller than 10~^, or when V is 
composed by less than 3000 particles. If the latter condition 
is not fulfilled the shape of such a contour is not considered 
in our analysis. 

A2 Density 

An alternative approach to determine halo shapes is to con- 
sider the underlying density field, which carries more infor- 
mation about the internal mass distribution of the haloes 
than the inertia tensor. We explore here two different imple- 
mentations. 



Jing & Suto (2002| determine the shape by fitting el- 



lipsoids to sets of particles having (nearly) the same nearest 
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Figure Al. Dark matter density in thin slices centred on the z = (left) y = (centre) and x = (right) planes. The black lines 
represent the isodensity contours. The solid white ellipses correspond to our new density-based method, while the dashed ones are for 
the inertia tensor method described by |Allgood et al. (2006). In the inner regions both methods agree very well and follow the isodensity 
surface, however closer to the virial radius (the external most contour) the inertia tensor-based scheme predicts somewhat rounder shapes 
than the actual density distribution. 



neighbours-based density. Noise due to substructures is ef- 
fectively removed by an implementation of the Friends of 



Friends (fof) algorithm to these sets ( [Davis et al. 1985), 
and where the linking length is selected to vary from set 
to set of iso-density particles according to the empirical law 
/ = 3(y9/m)-^/^ 

• We present a new method based on fitting an ellipsoid 
to the particle density. We first create a continuous den- 
sity field out of the particles positions. We do so by using 
a Cloud-In-Cell algorithm that allows the reconstruction of 
the density field on a regular grid ( Hockney Eastwood] 
1988). For better resolution and to keep the computational 



cost relatively low, the region covered by the grid is itera- 
tively increased. The second step involves the identification 
of iso-density contours. We select the cells with nearly the 
same density and a version of the fof algorithm is used to 
get rid of cells associated to substructures artificially linked 
to the main contour. Finally, we minimize the function: 



(A2) 



with Al the matrix representation of an ellipsoid, in order 
to determine the axis lengths (eigenvalues of M) and direc- 
tions (eigenvectors). Notice that the minimization is carried 
out in a 6 dimensional space (iVf has just 6 independent 
elements), therefore an educated initial guess for the itera- 
tion may reduce the numerical effort. We provide this guess 
by diagonalizing the inertia tensor of the cells with similar 
values of the density. 



Fig. |Al| shows 2D slices of the density map computed for 
halo Aq-A-4, together with several best-fit ellipsoids found 
by the method just outlined. An important feature of our 
method is that all the information about the 3D isoden- 
sity surface is taken into account; this is particularly useful 
towards the outskirts of the dark matter haloes where, as 
can be appreciated from Fig. |A1| density contours become 
less symmetric. For comparison, we also show the results 



of applying the algorithm by Allgood et al. (2006) (after 
subhaloes subtraction) in dashed lines. Recall that in our 
new method ellipsoids are effectively independent of each 



other and therefore this method is more sensitive to local 
variations of the halo shapes. On the other hand, |Allgood] 
et al. (2006) use the whole set of particles within a given 



radius, which implies that the shape of a contour at a given 
distance is correlated with the shape at smaller radii, as a 
careful inspection of this figure shows. 



A3 Potential 

A viable alternative to the density-based methods to mea- 
sure the shape of a dark matter halo is to use the grav- 
itational potential field defined by the particles. As first 



noticed by Springel et al. (2004) and later confirmed by 



Hayashi et al. (2007), substructures cause significant fiuc- 



tuations in the local density distribution, but their con- 
tribution to the gravitational potential is considerable less 



harmful (Springel et al. 



2004 



Hayashi et al. 2007). Iso- 



potential contours are therefore smoother and more regular 
than those defined by the density field. Taking advantage 



of this feature, Hayashi et al. (2007| have implemented a 
method to characterize the structural shape of dark mat- 
ter haloes by fitting ellipses to the iso-potential contours 
computed on three orthogonal 2D planes. It is important 
to note here that the potential is intrinsically rounder than 
the density distribution, for example for a cored logarithmic 
potential 1 — {c/a)p ~ 3[1 — (c/a)$] outside the core (Binney 
Tremaine||2008| ). 



A4 Results 



In Fig. |A2| we compare the shape as a function of the ellip- 
soidal radius obtained for the halo Aq-A-4 using the methods 
described above. This figure shows very good agreement in 
the halo shapes measured by the different methods, espe- 
cially in the inner regions {R < 100/i~^ kpc). There is how- 
ever an indication that at large radii, density-based methods 
(red solid and light-blue dashed curves) tend to give slightly 
more oblate shapes (higher b/a) than those based on im- 
plementations of the inertia tensor, but the effect is only 
marginal. As expected, the method based on the gravita- 
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Figure Bl. Distribution of particles around Aq-A halo at four 
different times. Particles with a single caustic crossing are plotted 
in red and are seen to trace reasonably well the filamentary struc- 
ture surrounding the halo. Each box has been rotated according 
to the inertia tensor defined by this subset of particles. 



tional potential ( [Hayashi et al.|[2007 ) produces higher axis 
ratios (blue dotted line). 



APPENDIX B: THE SIZE OF FILAMENTS AT 
DIFFERENT TIMES 

Dark matter haloes that reside in a filament will acquire 
their mass preferentially along the filament's longest axis. 
As discussed in Sec. [5] such infall of material gives rise to a 
power spectrum characterized by a significant / = 2 compo- 
nent. However, when the surrounding filament is sufficiently 
wide, i.e. of comparable or larger cross-section than the virial 
radius of the halo, the infalling particles will appear to be 
more isotropically distributed on the sky, shifting the power 
to the I — term of the spectrum. We have argued in Section 
|5]ofthis paper that the latter case is characteristic of the late 
stages of mass assembly in ^ lO^^M© objects. We analyze 
this statement in more detail in this Appendix, and provide 
a suitable measurement of a filament's size and compare it 
to that of the dark matter halo it hosts. For brevity we focus 
our analysis on halo Aq-A, but our conclusions should hold 
in the general case. 

The measurement of the size of a filament is not com- 
pletely straightforward and may depend on the particular al- 
gorithm used for its identification ( S toica et al.|2005 Zhang 
et al.||2009 ). In this work we define structures according to 



the number of caustic crossings of its constituent particles. 
Caustics arise during the gravitational collapse of a dynami- 
cal system. In the case of cold dark matter, the initial veloci- 
ties of particles are negligible, and thus these are distributed 
in 3D sheets in phase-space ( Bertschinger|| 19851 . Their col- 
lapse and subsequent virialization may be seen as the folding 
of these sheets, and the location of these folds in configura- 
tion space gives rise to caustics. Therefore the number of 



caustic crossings is indicative of the degree of virialization 
of a dynamical system. 

For example particles with zero crossings remain un- 
der the quasi-linear regime and have therefore not collapsed 
into any virialized structure today. As gravitational collapse 
proceeds, the number of caustics that a particle experiences 
increases rapidly. As shown in |Vogelsberger k, White| (2011 



Fig. 4), particles with 1 or 2 caustic crossings delineate the 
surrounding filamentary structure of a halo, whereas those 
with a higher number of crossings belong to the host halo 
itself. Therefore, in this work we select particles with 1 caus- 
tic crossing to study the properties of the filamentary struc- 



ture surrounding Aq-A halo at different epochs ( [Vogelsberger 
et al.||20Q9||Vogelsberger White||2011| ). 

Fig. |B1| shows the distribution of particles in a box of 
size 6rvir at four different times, where those with a single 
caustic crossing are highlighted in red. The circle indicates 
the corresponding virial radius of the Aq-A halo, and the hor- 
izontal bar gives the reference for conversion to the physical 
scale. At each time-step we have rotated the reference frame 
to the principal axis of the inertia tensor defined by the 
particles with 1 caustic crossing This figure shows these 
particles successfully trace the filamentary structure around 
this halo. Interestingly, there is a hint of evolution in the 
relative size of the filament with respect to that of the halo: 
as we move back in time, the red particles change from fully 
encompassing the halo (t ^ 9 Gyr) to having a similar cross- 
section (t ^ 5.25 Gyr) to becoming even narrower at earlier 
times (t ^ 2.3 Gyr). At this point the geometry of the sys- 
tem is much more complex due to multiple filaments feeding 
material to the central object. 

In order to quantify the evolution in the cross section 
of a filament we proceed as follows. We rotate the reference 
system such that the z-direction is that given by the iner- 
tia tensor of particles with 1 caustic crossing. We identify 
planes perpendicular to this direction and label the positions 
of such planes as Zp, where Zp — [—2, —1.5, —1, +1.5, +2] 
rvir- We project in each plane all selected particles that sat- 
isfy \zi — Zp\ ^ 0.5rvir- We compute the (projected) density of 
neighbours of each particle in the plane using a 2D SPH ker- 
nel and identify the particle with the highest density. This 
particle's location sets the centroid of the filament and its 
density at the core. We then define the radius of the filament 
rfii as the distance from this centroid where the density has 
dropped to 60% of its central value. Although this definition 
is arbitrary, this allows a measurement of the relative size 
at different times. 

Fig. [B2] illustrates our procedure for two different times: 
t = 13.6 (top) and t = 2.3 (bottom) Gyr; applied to 3 differ- 
ent planes along the z-direction: Zp — — rvir, 0, +rvir. Each 
panel shows the projected density of particles belonging to 
the filamentary structure (one single caustic crossing) , where 
the virial radius of the central halo is indicated with a white 
(empty) circle, and the highlighted full disc indicates the size 
of the filament centred on the highest density point as de- 



^ Because we are interested in the set of filaments connected to 
the central object, we run a FOF algorithm over the particles with 
1 caustic crossing, using a linking length of 0.7 the mean inter- 
particle separation and retain for the analysis only those set of 
particles belonging to the most massive FOF group 



16 Vera- Giro et al 



lOOfe"^ kpc 



13.58 Gyr ■ Orvir 



13.58 Gyr ■ l.rvir 





Figure B2. Projected mass density (in units of /iM0kpc~^) for three different perpendicular planes located along the direction of the 
filament at positions Zp = [—1,0, l]rvir, and at two different times: t = 13.6 Gyr (top) and t = 2.3 Gyr (bottom). The solid white curve 
indicates the half the virial radius of the halo at each time and the highlighted filled disc shows the size assigned to the filament and 
centred in the highest density point as described in the text. The geometry of the filaments is complicated, but their relative size with 
respect to that of the halo is clearly larger at present day {t ~ 13.6 Gyr) than it was at t ~ 2.3 Gyr. 



scribed above. Fig. |B2| suggests that filaments are not really 
straight in space (as indicated by the different positions of 
the centroids for panels located at different Zp) and also their 
geometry is complex since panels located symmetrically with 
respect to the centre of the halo yield significantly different 
sizes (e.g. top row, — rvir, +^vir panels). However, there is 
still a clear trend indicating that the filament's relative size 
was smaller at early times than at the present day. 

This trend is more clearly seen in Fig. |B3| where we 
show the average size of the filament (normalized to a 
fraction of the instantaneous virial radius) measured as 
the average over 9 equally distant planes located from 
[— 2rvir, +2rvir] as a function of time. The shaded region 
shows the scatter between the sizes derived for each of the 
individual planes at a fixed time. Notice that the relative 
filament- to- halo size at present day is almost twice as large 
than its value at t ^ 2 Gyr. This provides further support 
to our observation that the infall of particles at later times 
results from an increased size of the filament with respect 
to that of the halo. 




t (Gyr) 

Figure B3. Filament size relative to the virial radius / = 
rfii/(0.5rvir) as a function of time. The solid black line is the 
average size of the filament measured on 9 planes perpendicular 
to the direction of the filament. The shaded region is la scatter 
(see text for details) . The vertical axis has been normalized to its 
final value at redshift zero /o = f(z = 0) for an easier comparison. 



